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METHOD AND SYSTEM FOR SIMULATING A 



HYDROCARBON-BEARING FORMATION 



This application claims the benefit of U.S. Provisional Application No. 
60/159,035 filed on October 12, 1999. 



10 



FIELD OF THE INVENTION 



This invention relates generally to simulating a hydrocarbon-bearing 
formation, and more specifically to a method and system for simulating a 
hydrocarbon-bearing formation under conditions in which a fluid is injected into the 
formation to displace resident hydrocarbons. The method of this invention is 
1 5 especially useful in modeling the effects of viscous fingering and channeling as the 
injected fluid flows through a hydrocarbon-bearing formation. 



In the primary recovery of oil from a subterranean, oil-bearing formation or 
reservoir, it is usually possible to recover only a limited proportion of the original oil 

20 present in the reservoir. For this reason, a variety of supplemental recovery techniques 
have been used to improve the displacement of oil from the reservoir rock. These 
techniques can be generally classified as thermally based recovery methods (such as 
steam flooding operations), waterflooding methods, and gas-drive based methods that 
can be operated under either miscible or immiscible conditions. 

25 In miscible flooding operations, an injection fluid or solvent is injected into the 

reservoir to form a single-phase solution with the oil in place so that the oil can then 
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be removed as a more highly mobile phase from the reservoir. The solvent is typically 
a light hydrocarbon such as liquefied petroleum gas (LPG), a hydrocarbon gas 
containing relatively high concentrations of aliphatic hydrocarbons in the C2 to Ce 
range, nitrogen, or carbon dioxide. Miscible recovery operations are normally carried 
5 out by a displacement procedure in which the solvent is injected into the reservoir 
through an injection well to displace the oil from the reservoir towards a production 
well from which the oil is produced. This provides effective displacement of the oil in 
the areas through which the solvent flows. Unfortunately, the solvent often flows 
unevenly through the reservoir. 

10 Because the solvent injected into the reservoir is typically substantially less 

viscous than the resident oil, the solvent often fingers and channels through the 
reservoir, leaving parts of the reservoir unswept Added to this fingering is the 
inherent tendency of a highly mobile solvent to flow preferentially through the more 
permeable rock sections or to gravity override in the reservoir. 

1 5 The solvent's miscibility with the reservoir oil also affects its displacement 

efficiency within the reservoir. Some solvents, such as LPG, mix directly with 
reservoir oil in all proportions and the resulting mixtures remain single phase. Such 
solvent is said to be miscible on first contact or "first-contact miscible." Other 
solvents used for miscible flooding, such as carbon dioxide or hydrocarbon gas, form 

20 two phases when mixed directly with reservoir oil — therefore they are not first-contact 
miscible. However, at sufficiently high pressure, in-situ mass transfer of components 
between reservoir oil and solvent forms a displacing phase with a transition zone of 
fluid compositions that ranges from oil to solvent composition, and all compositions 
within the transition zone of this phase are contiguously miscible. Miscibility 

25 achieved by in-situ mass transfer of the components resulting from repeated contact of 
oil and solvent during the flow is called "multiple-contact" or dynamic miscibility. 
The pressure required to achieve multiple-contact miscibility is called the "minimum- 
miscibility pressure." Solvents just below the minimum miscibility pressure, called 
"near-mi scible" solvents, may displace oil nearly as well as miscible solvents. 

30 Predicting miscible flood performance in a reservoir requires a realistic model 

representative of the reservoir. Numerical simulation of reservoir models is widely 
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used by the petroleum industry as a method of using a computer to predict the effects 
of miscible displacement phenomena. In most cases, there is desire to model the 
transport processes occurring in the reservoir. What is being transported is typically 
mass, energy, momentum, or some combination thereof. By using numerical 
5 simulation, it is possible to reproduce and observe a physical phenomenon and to 
determine design parameters without actual laboratory experiments and field tests. 

Reservoir simulation infers the behavior of a real hydrocarbon-bearing 
reservoir from the performance of a numerical model of that reservoir. The objective 
is to understand the complex chemical, physical, and fluid flow processes occurring in 

10 the reservoir sufficiently well to predict future behavior of the reservoir to maximize 
hydrocarbon recovery. Reservoir simulation often refers to the hydrodynamics of 
flow within a reservoir, but in a larger sense reservoir simulation can also refer to the 
total petroleum system which includes the reservoir, injection wells, production wells, 
surface flowlines, and surface processing facilities. 

1 5 The principle of numerical simulation is to numerically solve equations 

describing a physical phenomenon by a computer. Such equations are generally 
ordinary differential equations and partial differential equations. These equations are 
typically solved using numerical methods such as the finite element method, the finite 
difference method, the finite volume method, and the like. In each of these methods, 

20 the physical system to be modeled is divided into smaller gridcells or blocks (a set of 
which is called a grid or mesh), and the state variables continuously changing in each 
gridcell are represented by sets of values for each gridcell. In the finite difference 
method, an original differential equation is replaced by a set of algebraic equations to 
express the fundamental principles of conservation of mass, energy, and/or 

25 momentum within each gridcell and transfer of mass, energy, and/or momentum 
transfer between gridcells. These equations can number in the millions. Such 
replacement of continuously changing values by a finite number of values for each 
gridcell is called "discretization". In order to analyze a phenomenon changing in time, 
it is necessary to calculate physical quantities at discrete intervals of time called 

30 timesteps, irrespective of the continuously changing conditions as a function of time. 



I:\URC\A-LAW\PATENTS\RES\97095VPROVISIONAL APPLICATION.DOC 




-4- 

Time-dependent modeling of the transport processes proceeds in a sequence of 
timesteps. 

In a typical simulation of a reservoir, solution of the primary unknowns, 
typically pressure, phase saturations, and compositions, are sought at specific points in 
5 the domain of interest. Such points are called "gridnodes 11 or more commonly 

"nodes." Gridcells are constructed around such nodes, and a grid is defined as a group 
of such gridcells. The properties such as porosity and permeability are assumed to be 
constant inside a gridcell. Other variables such as pressure and phase saturations are 
specified at the nodes. A link between two nodes is called a "connection." Fluid flow 

10 between two nodes is typically modeled as flow along the connection between them. 

Compositional modeling of hydrocarbon-bearing reservoirs is necessary for 
predicting processes such as first-contact miscible, multiple-contact miscible, and 
near-miscible gas injection. The oil and gas phases are represented by 
multicomponent mixtures. In such modeling, reservoir heterogeneity and viscous 

15 fingering and channeling cause variations in phase saturations and compositions to 
occur on scales as small as a few centimeters or less. A fine-scale model can 
represent the details of these adverse-mobility displacement behaviors. However, use 
of fine-scale models to simulate these variations is generally not practical because 
their fine level of detail places prohibitive demands on computational resources. 

20 Therefore, a coarse-scale model having far fewer gridcells is typically developed for 
reservoir simulation. Considerable research has been directed to developing models 
suitable for use in predicting miscible flood performance. 

Development of a coarse-grid model that effectively simulates gas 
displacement processes is especially challenging. For compositional simulations, the 

25 upscaled, coarse-grid model must effectively characterize changes in phase behavior 
and changes in oil and gas compositions as the oil displacement proceeds. Many 
different techniques have been proposed. Most of these proposals use empirical 
models to represent viscous fingering in first-contact miscible displacement. See for 
example: 
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Koval, E. J., "A Method for Predicting the Performance of Unstable Miscible 
Displacement in Heterogeneous Media, " Society of Petroleum Engineering 
Journal, pages 145-154, June 1963; 

Dougherty, E. L., "Mathematical Model of an Unstable Miscible Displacement," 
5 Society of Petroleum Engineering Journal, pages 155-163, June 1963; 

Todd, M. R., and Longstaff, W. J., "The Development, Testing, and Application of 
a Numerical Simulator for Predicting Miscible Flood Performance," Journal of 
m Petroleum Technology, pages 874-882, July 1972; 

%Q Fayers, F. J., "An Approximate Model with Physically Interpretable Parameters 

1 0 for Representing Miscible Viscous Fingering," SPE Reservoir Engineering, pages 

l 5 542-550, May 1988; and 

C3 Fayers, F. J. and Newley, T. M. J., "Detailed Validation of an Empirical Model for 

" Viscous Fingering with Gravity Effects," SPE Reservoir Engineering, pages 542- 

O 550, May 1988. 

fO 15 Of these models, the Todd-Longstaff ("T-L") mixing model is the most 

VZ popular, and it is used widely in reservoir simulators. When properly used, the T-L 

3 mixing model provides reasonably accurate average characteristics of adverse- 

mobility displacements when the injected solvent and oil are first-contact miscible. 
However, the T-L mixing model is less accurate under multiple-contact miscible 
20 conditions. 

Models have been suggested that use the T-L model to account for viscous 
fingering under multiple-contact miscible situations (see for example Todd, M. R. and 
Chase, C. A., "A Numerical Simulator for Predicting Chemical Flood Performance," 
SPE-7689, presented at the 54th Annual Fall Technical Conference and Exhibition of 

25 the Society of Petroleum Engineers, Houston, Texas, 1979, sometimes referred to as 
the "Todd-Chase technique"). In modeling a multiple-contact miscible displacement, 
in addition to the viscous fingering taken into account in the T-L mixing model, 
exchange of solvent and oil components between phases according to the phase 
behavior relations must also be considered. The importance of the interaction between 

30 phase behavior and fingering in multiple-contact miscible displacements was disclosed 
by Gardner, J. W., and Ypma, J. G. J., "An Investigation of Phase- 
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Behavior/Macroscopic Bypassing Interaction in C0 2 Flooding, " Society of Petroleum 
Engineering Journal, pages 508-520, October 1984. However, these proposals did not 
effectively combine use of a mixing model and a phase behavior model. 

Another proposed model for taking into account fingering and channeling 
5 behavior in multiple-contact miscible displacement suggested making the 

dispersivities of solvent and oil components dependent on the viscosity gradient, 
thereby addressing the macroscopic effects of viscous fingering (see Young, L. C, 
"The Use of Dispersion Relationships to Model Adverse Mobility Ratio Miscible 
Displacements," paper SPE/DOE 14899 presented at the 1986 SPE/DOE Enhanced 

10 Oil Recovery Symposium, Tulsa, April 20-23). Another model proposed extending 
the T-L model to multiphase multicomponent flow with simplified phase behavior 
predictions (see Crump, J. G., "Detailed Simulations of the Effects of Process 
Parameters on Adverse Mobility Ratio Displacements," paper SPE/DOE 17337, 
presented at the 1988 SPE/DOE Enhanced Oil Recovery Symposium, Tulsa, 

15 April 17-20). A still another model suggested using the fluid compositions flowing 
out of a large gridcell to compensate for the nonuniformity of the fluid distribution in 
the gridcell (see Barker, J. W., and Fayers, F. J., "Transport Coefficients for 
Compositional Simulation with Coarse Grids in Heterogeneous Media", SPE 22591, 
presented at SPE 66th Annual Tech. Conf., Dallas, TX, Oct. 6-9, 1991). A still 

20 another model proposed that incomplete mixing between solvent and oil can be 

represented by assuming that thermodynamic equilibrium prevails only at the interface 
between the two phases, and a diffusion process drives the oil and solvent composition 
towards these equilibrium values (see Nghiem, L. X., and Sammon, P. H., "A Non- 
Equilibrium Equation-of- State Compositional Simulator," SPE 37980, presented at the 

25 1997 SPE Reservoir Simulation Symposium, Dallas, TX, June 8-17, 1997). The 
gridcells in these models were not subdivided. 

Proposals have been made to represent fingering and channeling in multiple- 
contact miscible displacements using two-region models. See for example: 

Nghiem, L. X., Li, Y. K. and Agarwal, R. K., "A Method for Modeling Incomplete 

30 Mixing in Compositional Simulation of Unstable Displacements," SPE 18439, 
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presented at the 1989 Reservoir Simulation Symposium, Houston, TX, February 6- 
8, 1989; and 

Fayers, F. J., Barker, J. W., and Newley, T. M. J., "Effects of Heterogeneities on 
Phase Behavior in Enhanced Oil Recovery," in The Mathematics of Oil Recovery, 
5 P. R. King, editor, pages 1 15-150, Clarendon Press, Oxford, 1992. 

These models divide a simulation gridcell into a region where complete mixing occurs 
between the injected solvent and a portion of the resident oil and a region where the 
resident oil is bypassed and not contacted by the solvent. Although the conceptual 
structure of these models appears to provide a better representation of incomplete 

1 0 mixing in multiple-contact miscible displacements than single zone models, the 

physical basis of the equations used to represent bypassing and mixing is unclear. In 
particular, these models (1) use empirical correlations to represent oil/solvent 
mobilities in each region, (2) use empirical correlations to represent component 
transfer between regions, and (3) make restrictive assumptions about the composition 

15 of the regions and direction of component transfer between the regions. It has been 
suggested that the empirical mobility and mass transfer functions in these models can 
be determined by fitting them to the results of fine-grid simulations. As a result, in 
practice, calibration of these models will be a time-consuming and expensive process. 
Furthermore, these models are unlikely to accurately predict performance outside the 

20 parameter ranges explored in the reference fine-grid simulations. 

While the two-region approaches proposed in the past have certain advantages, 
there is a continuing need for improved simulation models that provide a better 
physical representation of bypassing and mixing in adverse mobility displacement and 
thus enable more accurate and efficient prediction of flood performance. 

25 

SUMMARY 

A method and system is provided for simulating one or more characteristics of 
a multi-component, hydrocarbon-bearing formation into which a displacement fluid 
having at least one component is injected to displace formation hydrocarbons. The 
30 first step of the method is to equate at least part of the formation to a multiplicity of 
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gridcells. Each gridcell is then divided into two regions, a first region representing a 
portion of each gridcell swept by the displacement fluid and a second region 
representing a portion of each gridcell essentially unswept by the displacement fluid. 
The distribution of components in each region is assumed to be essentially uniform. A 
5 model is constructed that is representative of fluid properties within each region, fluid 
flow between gridcells using principles of percolation theory, and component transport 
between the regions. The model is then used in a simulator to simulate one or more 
characteristics of the formation. 

1 0 BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention and its advantages will be better understood by referring 
to the following detailed description and the following drawings in which like 
numerals have similar functions. 

Fig. 1 illustrates a two-dimensional schematic of a solvent flowing through an 
1 5 oil reservoir to displace oil therefrom, which shows an example of solvent fingering 
in the reservoir. 

Fig. 2 illustrates an example of a two-dimensional fine-scale grid that could 
represent the reservoir area of Fig. 1. 

Fig. 3 illustrates a two-dimensional gridcell covering the same domain 
20 depicted in Fig. 1, with the gridcell divided into two regions, one representing the 
region of the domain swept by an injected fluid and the second region representing 
the region of the domain unswept by the injected fluid. 

Fig. 4 illustrates the gridcell depicted in Fig. 3 showing schematically phase 
fractions in the two regions of the gridcell. 
25 Fig. 5 A illustrates the effect of coordination number, z, on total oil recovery 

for a multiple-contact miscible flood simulated using the method of this invention. 

Fig. 5B illustrates the effect of coordination number, z, on solvent 
breakthrough for a multiple-contact miscible flood simulated using the method of this 
invention. 
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Figs. 6A-D illustrate the effect of oil Damkohler numbers on heavy and light 
oil recovery curves for a multiple-contact miscible flood simulated using the method 
of this invention. 

Fig. 7 graphically compares published first-contact miscible flood recovery 
5 data and best-fits obtained using the method of this invention. 

Fig. 8 illustrates coordination numbers obtained by fitting the model used in 
the method of this invention and published data as a function of oil/solvent viscosity 
ratio. 

Fig. 9 illustrates published experimental C0 2 /Soltrol and C0 2 /Wasson crude 
1 0 coreflood recovery data and simulation predictions using a published single-region 
model. 

Fig. 10 illustrates published experimental CCVSoltrol and CCVWasson crude 
coreflood recovery data and simulation predictions using the method of this invention. 

The drawings illustrate specific embodiments of practicing the method of this 
15 invention. The drawings are not intended to exclude from the scope of the invention 
other embodiments that are the result of normal and expected modifications of the 
specific embodiments. 



DETAILED DESCRIPTION OF THE INVENTION 

In order to more fully understand the present invention, the following 
20 introductory comments are provided. To increase the recovery of hydrocarbons from 
subterranean formation, a variety of enhanced hydrocarbon recovery techniques have 
been developed whereby a fluid is injected into a subterranean formation at one or 
more injection wells within a field and hydrocarbons (as well as the injected fluid) are 
recovered from the formation at one or more production wells within the field. The 
25 injection wells are typically spaced apart from the production wells, but one or more 
injection wells could later be used as production wells. The injected fluid can for 
example be a heating agent used in a thermal recovery process (such as steam), any 
essentially immiscible fluid used in an immiscible flooding process (such as natural 
gas, water, or brine), and any miscible fluid used in a miscible flooding process (for 
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example, a first-contact miscible fluid, such as liquefied petroleum gas, or a multiple- 
contact miscible or near-miscible fluid such as lower molecular weight hydrocarbons, 
carbon dioxide, or nitrogen). 

Fig. 1 schematically illustrates a two-dimensional reservoir area 5 which is part 
5 of a larger oil-bearing, geologic formation (not shown) to be analyzed using the 
method of this invention. In Fig. 1, an injected fluid 1 1, which is assumed to be 
gaseous in this description, displaces a multi-component resident oil 12 in the 
reservoir area 5. It should be understood that this invention is not limited to a gaseous 
injected fluid; the injected fluid could also be liquid or a multi-phase mixture. The 

10 injected fluid 1 1 flows from left to right in the drawing. Fig. 1 depicts viscous 

fingering that occurs when the injected fluid 1 1 displaces resident oil 12. The injected 
fluid 1 1 tends to finger through the oil 12 towards a producing well (not shown in the 
drawing), resulting in premature breakthrough of the injected fluid 1 1 and bypassing 
some of the resident oil 12. Viscous fingering is dominantly caused by large 

15 differences in oil 12 and injected fluid 1 1 viscosities resulting in a mobility ratio of 
injected fluid-to-oil that has an adverse effect on areal sweep efficiency or 
displacement efficiency of the injected fluid. 

Through advanced reservoir characterization techniques, the reservoir area 5 
can be represented by gridcells on a scale from centimeters to several meters, 

20 sometimes called a fine-scale grid. Each gridcell can be populated with a reservoir 
property, including for example rock type, porosity, permeability, initial interstitial 
fluid saturation, and relative permeability and capillary pressure functions. 

Fig. 2 shows an example of a two-dimensional fine-scale grid 10 that could 
represent the reservoir area 5 of Fig. 1. The reservoir area 5 of Fig. 1 is represented in 

25 Fig. 2 by 84 gridcells. Gridcells IT represent the geologic regions that have been 
swept by injected fluid 1 1 and the gridcells 12' represent the geologic regions that 
contain essentially resident oil 12 undisplaced by the injected fluid. However, 
reservoir simulations are not typically performed with such fine-scale grids. The 
direct use of fine-scale models for full-field reservoir simulation is not generally 

30 feasible because their fine level of detail places prohibitive demands on computational 
resources. Therefore, a coarse-scale grid is typically used in simulation models, while 
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preserving, as much as possible, the fluid flow characteristics and phase behavior of 
the fine-scale grid. A coarse-scale grid may represent, for example, all 84 gridcells of 
Fig. 2 by one gridcell. A method is therefore needed to model fluid compositions and 
fluid flow behavior taking into account fingering and channeling. The method of this 
5 invention provides this capability. 

The method of this invention begins by equating the reservoir area to be 
analyzed to a suitable grid system. The reservoir area to be analyzed is represented by 
a multiplicity of gridcells, arranged adjacent to one another so as to have a boundary 
between each pair of neighboring gridcells. This spatial discretization of the reservoir 

1 0 area can be performed using finite difference, finite volume, finite element, or similar 
well-known methods that are based on dividing the physical system to be modeled into 
smaller units. The present invention is described primarily with respect to use of the 
finite difference method. Those skilled in the art will recognize that the present 
invention can also be applied in connection with finite element methods or finite 

15 volume methods. When using the finite difference and finite volume methods, the 
smaller units are typically called gridcells, and when using the finite element method 
the smaller units are typically called elements. These gridcells or elements can 
number from fewer than a hundred to millions. In this patent, for simplicity of 
presentation, the term gridcell is used, but it should be understood that if a simulation 

20 uses the finite element method the term element would replace the term gridcell as 
used in this description. 

In the practice of this invention, the gridcells can be of any geometric shape, 
such as parallelepipeds (or cubes) or hexahedrons (having four vertical corner edges 
which may vary in length), or tetrahedra, rhomboids, trapezoids, or triangles. The grid 

25 can comprise rectangular gridcells organized in a regular, structured pattern (as 

illustrated in Fig. 2), or it can comprise gridcells having a variety of shapes laid out in 
an irregular, unstructured pattern, or it can comprise a plurality of both structured and 
unstructured patterns. Completely unstructured grids can be assembled that assume 
almost any shape. All the gridcells are preferably boundary aligned, thereby avoiding 

30 having any side of a gridcell contacting the sides of two other gridcells. 
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One type of flexible grid that can be used in the model of this invention is the 
Voronoi grid. A Voronoi gridcell is defined as the region of space that is closer to its 
node than to any other node, and a Voronoi grid is made of such gridcells. Each 
gridcell is associated with a node and a series of neighboring gridcells. The Voronoi 
5 grid is locally orthogonal in a geometrical sense; that is, the gridcell boundaries are 
normal to lines joining the nodes on the two sides of each boundary. For this reason, 
Voronoi grids are also called perpendicular bisection (PEBI) grids. A rectangular grid 
block (Cartesian grid) is a special case of the Voronoi grid. The PEBI grid has the 
flexibility to represent widely varying reservoir geometry, because the location of 

10 nodes can be chosen freely. PEBI grids are generated by assigning node locations in a 
given domain and then generating gridcell boundaries in a way such that each gridcell 
contains all the points that are closer to its node location than to any other node 
location. Since the inter-node connections in a PEBI grid are perpendicularly bisected 
by the gridcell boundaries, this simplifies the solution of flow equations significantly. 

15 For a more detailed description of PEBI grid generation, see Palagi, C.L. and Aziz, 

K.: "Use of Voronoi Grid in Reservoir Simulation," paper SPE 22889 presented at the 
66th Annual Technical Conference and Exhibition, Dallas, TX (Oct. 6-9, 1991). 

The next step in the method of this invention is to divide each gridcell that has 
been invaded by the injected fluid into two regions, a first region that represents a 

20 portion of the gridcell swept by the injected fluid 1 1 and a second region that 
represents a portion of the gridcell that is unswept by the injected fluid 1 1 . The 
distribution of components in each region is assumed to be uniform. It is further 
assumed that fluids within each region are at thermodynamic equilibrium. However, 
the two regions of the gridcell are not in equilibrium with each other, and as a result 

25 the compositions and phase volume fractions within each region will typically be 
different. 

Fig. 3 illustrates a two-dimensional schematic of one grid gridcell 1 5 that 
represents the same reservoir area represented by the 84 gridcells of grid 10 (Fig. 2). 
While not shown in the drawing, it should be understood that gridcell 1 5 shares 
30 boundaries with neighboring gridcells. The following description with respect to 
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gridcell 15 also applies to other gridcells in the grid of which gridcell 15 is only one of 
a multiplicity of gridcells. 

Referring to Fig. 3, gridcell 15 is divided into two regions 16 and 17. 
Region 16 represents the portion of the gridcell invaded by the injected fluid 1 1 and 
5 region 1 7 represents the portion of the gridcell that has not been displaced by the 

injected fluid 11. Regions 16 and 17 are separated by an interface or partition 18 that 
is assumed to have infinitesimal thickness. Multicomponent fluids within each region 
are assumed to be in thermodynamic equilibrium, which means that the fluid 
compositions and phase volumes of regions 16 and 17 could be different, and typically 
; 10 are different. The compositions of fluids can vary from gridcell to gridcell within the 
in grid and the compositions of fluids within each region of a gridcell can vary with time. 

?5 Therefore, partition 18 can move as a function of time as the injected fluid 1 1 contacts 

^ more of the region represented by gridcell 15. Movement of partition 18 depends 

3 primarily on (1) exchange of fluids between gridcell 15 and its neighboring gridcells, 

;H 1 5 (2) mass transfer across the partition 1 8, and (3) injection or withdrawal of fluids 
□ through injection and production wells that may penetrate the geologic region 

3 represented by gridcell 15. 

Fig. 4 illustrates an example of phase fractions of fluids in regions 16 and 17. 
The fraction of vapor phase, which consists of the injected fluid plus vaporized oil, is 
20 shown by numeral 1 la in region 16 and by numeral 1 lb in region 17. The fraction of 
liquid phase, which consists of resident oil plus dissolved injected fluid, is shown by 
numeral 12a in region 16 and by numeral 12b in region 17. The fraction of water is 
shown by numeral 13a in region 16 and numeral 13b in region 17. In the example 
shown in Fig. 4, region 16 contains primarily the high-mobility injected fluid 1 1 and 
25 region 17 contains primarily the low-mobility resident oil 12. Arrow 20 represents a 
fluid stream flowing into region 16 from invaded regions of gridcells adjacent to 
gridcell 15. Arrow 21 represents a fluid stream flowing into region 17 from resident 
regions of gridcells adjacent to gridcell 15. Arrow 22 represents a fluid stream 
flowing out of region 16 into invaded regions of gridcells adjacent to gridcell 15. 
30 Arrow 23 represents a fluid stream flowing out of region 17 into resident regions of 
gridcells adjacent to gridcell 15. Although the arrows show fluid flowing from left to 
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right, the fluid could flow into and out of gridcell 15 in other directions. Arrows 24 
represent mass transfer between regions 16 and 17. Components are allowed to 
transfer in either direction across the partition 18. Although the arrows 24 show 
transfer between phases of the same type (vapor to vapor, liquid hydrocarbon to liquid 
5 hydrocarbon, and water to water), components may transfer from any phase in the 
source region into any phase in the other region. Region 16 has zero volume until 
injected fluid 1 1 flows into gridcell 15. Injected fluid 1 1 may be modeled as being 
injected into either the invaded region 16 or resident region 17, or the injected fluid 1 1 
may be modeled as being injected simultaneously into both regions 16 and 17. Fluids 

10 may be withdrawn from both invaded region 16 and resident region 17. Gridcell 15 
can also be modeled as having injected fluid 1 1 flowing from one or more injection 
wells directly into gridcell 15, and it can be modeled as having fluid flowing directly 
out of gridcell 15 into one or more production wells. Although not shown in the 
drawings, if the reservoir area represented by gridcell 15 is penetrated by an injection 

15 well, injected fluid 1 1 injected into gridcell 15 could be modeled as being injected 
only into the invaded region 16 and if the reservoir area represented by gridcell 15 is 
penetrated by a production well, gridcell 1 5 could be modeled as having fluids being 
produced from both invaded region 16 and resident region 17. 

Although the drawings do not show gridcell nodes, persons skilled in the art 

20 would understand that each gridcell would have a node. In simulation operations, 
flow of fluid between gridcells would be assumed to take place between gridcell 
nodes, or, stated another way, through inter-node connections. In practicing the 
method of this invention, the invaded region of a given gridcell (region 16 of Figs. 3 
and 4) is connected to invaded regions of gridcells adjacent to the given gridcell, and 

25 the resident regions of a given gridcell (region 17 of Fig. 2) is connected to resident 

regions of gridcells adjacent to the given gridcell. There are no inter-node connections 
between resident region 16 and invaded region 17. The inventors therefore sometimes 
refer to the method of this invention as the Partitioned Node Model or PNM. 

The next step in the method of this invention is to construct a predictive model 

30 that represents fluid properties within each region of each gridcell, fluid flow between 
each gridcell and its neighboring gridcells, and component transport between regions 
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16 and 17 for each gridcell. In a preferred embodiment, the model comprises a set of 
finite difference equations for each gridcell having functions representative of the 
mobility of each fluid phase in regions 16 and 17, functions representative of the phase 
behavior within regions 16 and 17, and functions representative of the mass transfer of 
5 each component between the regions 16 and 17. The model may optionally further 
contain functions representing energy transfer between regions 16 and 17. Energy 
transfer functions may be desired for example to simulate the heat effects resulting 
from a steam flooding operation. 

Mobility functions are used to describe flow through the connections, and a 

10 mobility function is generated for each phase in each region. The mobilities of the 
streams 22 and 23 leaving the gridcell 1 5 depend on many factors including the 
composition of the fluids in the invaded region 16 and the resident region 17, the 
relative sizes (or volume fraction) of the invaded region 16 and resident region 17, the 
heterogeneity of the gridcell, and the oil-to-injected fluid mobility ratio. The specific 

15 functional dependencies are determined through the use of percolation theory. The 

basic principles of percolation theory are described by S. Kirkpatrick, "Percolation and 
Conduction," Rev. Modern Physics, vol. 45, pages 574-588, 1973, which is 
incorporated herein by reference. In a preferred embodiment, an effective medium 
mobility model represents the gridcell by a pore network so as to characterize the 

20 effect of fingering and channeling that occurs in the gridcell depending on conditions 
prevalent in the gridcell over a time interval. The effective mobility of each fluid 
phase in each region of a gridcell can be calculated by those skilled in the art having 
benefit of the teaching of this description. Examples of phase mobility equations, 
derived from an effective medium model, are provided below as equations (18) - (20). 

25 The method of this invention assumes that equilibrium exists within the 

invaded region 16 and within resident region 17. As part of the model, a 
determination is made of the properties of the phases that coexist within regions 1 6 
and 17. Preferably, a suitable equation of state is used to calculate the phase behavior 
of region 16 and region 17. In the examples provided below, a one-dimensional 

30 model uses a simplified pseudoternary phase behavior model that characterizes 

mixtures of solvent and oil in terms of three pseudocomponents, solvent (CO2X a light 
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oil component, and a heavy oil component. The simplified phase behavior model is 
capable of simulating the salient features of displacements involving different degrees 
of miscibility ranging from first contact miscible, through multiple-contact miscible, 
and near-miscible, to immiscible. The phase behavior properties can be determined by 
5 persons of ordinary skill in the art. 

The method of this invention does not assume equilibrium between the invaded 
region 16 and the resident region 17 of a gridcell. Mass transfer functions are used to 
describe the rate of movement of components across the interface or partition 1 8 
between regions 16 and 17. This mass transfer is depicted in Fig. 4 by arrows 24. 

10 Mechanisms of mass transfer include, but are not limited to, molecular diffusion, 

convective dispersion, and capillary dispersion. The method of this invention assumes 
that each component's rate of mass transfer is proportional to a driving force times a 
resistance. Examples of driving forces include, but are not limited to, composition 
differences and capillary pressure differences between the two regions. Once a mass 

1 5 transfer function is generated for each fluid component, the rates of mass transfer 
depend on factors, including, but not limited to, component identity, degree of 
miscibility between the gas and oil, size of each region, gridcell geometry, gas/oil 
mobility ratio, velocity, heterogeneity, and water saturation. These functionalities can 
be built into the mass transfer model by those skilled in the art. Examples of mass 

20 transfer functions are provided as equations (10) and (14) - (16) below. 

One of the first steps in designing the model is to select the number of space 
dimensions desired to represent the geometry of the reservoir. Both external and 
internal geometries must be considered. External geometries include the reservoir or 
aquifer limits (or an element of symmetry) and the top and bottom of the reservoir or 

25 aquifer (including faults). Internal geometries comprises the areal and vertical extent 
of individual permeability units and non-pay zones that are important to the solution of 
the problem and the definition of well geometry (for example, well diameter, 
completion interval, and presence of hydraulic fractures emanating from the well). 
The model of this invention is not limited to a particular number of 

30 dimensions. The predictive model can be constructed for one-dimensional (1-D), two- 
dimensional (2-D), and three-dimensional (3-D) simulation of a reservoir. A 1-D 
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model would seldom be used for reservoir-wide studies because it can not model areal 
and vertical sweep. A 1-D gas injection model to predict displacement efficiencies 
can not effectively represent gravity effects perpendicular to the direction of flow. 
However, 1-D gas injection models can be used to investigate the sensitivity of 
reservoir performance to variations in process parameters and to interpret laboratory 
displacement tests. 

2- D areal fluid injection models can be used when areal flow patterns dominate 
reservoir performance. For example, areal models normally would be used to 
compare possible well patterns or to evaluate the influence of areal heterogeneity on 
reservoir behavior. 2-D cross-sectional and radial gas injection models can be used 
when flow patterns in vertical cross-sections dominate reservoir performance. For 
example, cross-sectional or radial models normally would be used to model gravity 
dominated processes, such as crestal gas injection or gas injection into reservoirs 
having high vertical permeability, and to evaluate the influence of vertical 
heterogeneity on reservoir behavior. 

3- D models may be desirable to effective represent complex reservoir 
geometry or complex fluid mechanics in the reservoir. The model can for example be 
a 3-D model comprising layers of PEBI grids, which is sometimes referred to in the 
petroleum industry as 2 1 / 2 -D. The layered PEBI grids are unstructured areally and 
structured (layered) vertically. Construction of layered 3-D grids is described by 

(1) Heinemann, Z. E., et al., "Modeling Reservoir Geometry With Irregular Grids," 
SPE Reservoir Engineering, May, 1991 and (2) Verma, S., et al., "A Control Volume 
Scheme for Flexible Grids in Reservoir Simulation," SPE 37999, SPE Reservoir 
Simulation Symposium, Dallas, TX, June, 1997. 

The present invention is not limited to dividing a gridcell into only two zones. 
The method of this invention could be used with gridcells having multiple partitions, 
thus dividing the gridcells into three or more zones. For example, a three-zone 
gridcell may have one zone representing the region of the reservoir invaded by an 
injected fluid, a second zone representing the region of the reservoir uninvaded by the 
injected fluid, and a third zone representing a mixing region of the reservoir's resident 
fluid and the injected fluid. In another example, in a steam injection operation, one 
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zone may represent the region of the reservoir invaded by the injected steam, a second 
zone may represent the region of the reservoir occupied by gas other than steam, and a 
third zone may represent the region of the reservoir not occupied by the injected steam 
or the other gas. The gas other than steam could be, for example, solution gas that has 
5 evolved from the resident oil when the reservoir pressure falls below the bubble point 
of the oil, or a second injected gas such as enriched gas, light hydrocarbon gas, or 



The method of this invention can be used to simulate recovery of oil from 
viscous oil reservoirs in which thermal energy is introduced into the reservoir to heat 

10 the oil, thereby reducing its viscosity to a point that the oil can be made to flow. The 
thermal energy can be in a variety of forms, including hot waterflooding and steam 
injection. The injection can be in one or more injection wells and production of oil 
can be through one or more spaced-apart production wells. One well can also be used 
for both injection of fluid and production of oil. For example, in the "huff and puff 

1 5 process, steam is introduced through a well (which can be a vertical or horizontal well) 
into a viscous hydrocarbon deposit for a period of time, the well is shut in to permit 
the steam to heat the hydrocarbon, and subsequently the well is placed on production. 

Once the predictive model is generated, it can be used in a simulator to 
simulate one or more characteristics of the formation as a function of time. The basic 

20 flow model consists of the equations that govern the unsteady flow of fluids in the 

reservoir grid network, wells, and surface facilities. Appropriate numerical algorithms 
can be selected by those skilled in the art to solve the basic flow equations. Examples 
of numerical algorithms that can be used are described in Reservoir Simulation, Henry 
L. Doherty Series Monograph, Vol. 13, Mattax, C. C. and Dalton, R. L., editors, 

25 Society of Petroleum Engineers, Richardson, TX, 1990. The simulator is a collection 
of computer programs that implement the numerical algorithms on a computer. 

Persons skilled in the art will readily understand that the practice of the present 
invention is computationally intense. Accordingly, use of a computer, preferably a 
digital computer, to practice the invention is virtually a necessity. Computer programs 

30 for various portions of the modeling process are commercially available (for example, 
software is commercially available to develop gridcells, display results, calculate fluid 
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flow properties, and solve linear set of equations that are used in a simulator). 
Computer programs for other portions of the invention could be developed by persons 
skilled in the art based on the teachings set forth herein. 

The practice of this invention can be applied to part or all gridcells in a grid 
5 system being modeled. To economize on computational time, the additional 

computations associated with dividing gridcells into two or more zones is preferably 
applied only to those gridcells simulation model that are being invaded by injected 
fluid. 

The method of this invention is an improvement over two-region displacement 
10 models used in the past. This improvement can be attributed to the following key 

differences. First, percolation theory is used to characterize the effect of fingering and 
channeling on effective fluid mobilities. Second, the rate of component transfer 
between regions is proportional to a driving force times a resistance. Third, the mass 
transfer functions account for actual mixing processes such as molecular diffusion, 
1 5 convective dispersion, and capillary dispersion. These improvements result in more 
accurate and efficient prediction of adverse mobility displacements. 

One-Dimensional Simulation Examples 

A one-dimensional model of this invention was generated and the model was 
20 tested using a proprietary simulator. Commercially available simulators could be 
readily modified by those skilled in the art using the teachings of this invention and 
the assumptions presented herein to produce substantially similar results to those 
presented below. In the model, allocation of components between resident and 
invaded regions was determined by transport equations that accounted for convection 
25 of the invaded and resident fluids and the rate of each component's transfer between 
the regions. A four-component fluid description was used in the simulator. The four 
components were solvent (CO2X a light fraction of crude oil, a heavy fraction of crude 
oil, and water. It was assumed that the fluids were incompressible and that ideal 
mixing occurred, which allowed the pressure equations to be de-coupled from the 
30 component transport equations and substitution of volume fractions for mole fractions 
as the compositional variables. Persons skilled in the art would be familiar with 
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15 



techniques of accounting for fluid compressibilities and non-ideal mixing. It was also 
assumed that the solvent did not transfer into the resident region and that water 
saturation was the same in both regions. 

The following description of the simulation examples refers to equations 
having a large number of mathematical symbols, many of which are defined as they 
occur throughout the text. Additionally, for purposes of completeness, a table 
containing definitions of symbols used herein is presented following the detailed 
description. 

The simulator was formulated in terms of the standard transport equations for 
the total amount of each component, augmented by transport equations for the amount 
of each component in the resident region. The amount of each component in the 
invaded region was then obtained by difference. Under these assumptions, the 
dimensionless transport equations for total solvent, heavy component of the oil, and 
water were, respectively: 



dw x 



d_ 



ile X \ + ^roe X r\ 



l) 



(l) 



dw 2 



d_ 



ite X 2 + ^roe X r2 



-1) 



A, 



(2) 



20 



dx 



d_ 



' w as 



(3) 



The total light component volume fraction, W3, was obtained from: 



25 



w 3 = 1 - w x - w 2 - S„ 



(4) 



In Eq. (4), component 1 is solvent, component 2 is the heavy fraction of the oil, and 
component 3 is the light fraction of the oil. 
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In Eqs. (1) through (4), g^x/L, r= ut/<j>L, J3= k/uL, A t = A ive + X iie + A roe + A w , 
L is core length, k is permeability, 0 is porosity, /? c is the capillary pressure between oil 
and water, yj is the volume fraction of component j in the vapor portion of the invaded 
region, Xj is the volume fraction of component j in the liquid portion of the invaded 
region, and x rj is the volume fraction of component j in the nonaqueous portion of the 
resident region. Wj = w rJ + Wy is the total volume fraction of component j, where Wy = 
0(Sgyj + SiXj) is the volume fraction of component j in the invaded region and w rj = (1 
- 6)(1 -S w )x rJ is the volume fraction of component j in the resident region. 6 is the 
volume fraction of the invaded region, defined as: 



(5) 



Sg and Si are, respectively, the vapor and liquid saturations in the invaded region. Xr 0e 
is the mobility of the resident fluid, X ive is the mobility of the vapor phase in the 
invaded region, Au e is the mobility of the liquid phase in the invaded region, and A w is 
the mobility of water, all calculated using effective medium theory, as described 
below. The total injection velocity, u, was assumed to be constant. 

The dimensionless transport equations for resident solvent, heavy oil, and light 
oil were, respectively: 



chv, 



r\ 



dr 



at 



(6) 



dr 



d_ 

at 



^ roe X r2 (/^ w ~ 1) 

A, 



A 2 ^Z 



(7) 
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A ^ L (8) 



where Aj is the rate of transfer (volume/time) of component j from the resident to 
invaded region. The first term on the right side of these equations accounted for 
5 convection of each component within the resident region, and the second term 
accounted for transfer of each component from the resident region to the invaded 
region. 

The equation for pressure was: 

io '* ' ( 9) 

In the one-dimensional simulator, equations (1) through (3) and (6) through (8) were 
discretized to produce six sets of finite-difference equations in £ which are solved 
time-wise with Hamming's predictor-corrector method of integrating a set of first- 
order ordinary differential equations (the Hamming method would be familiar to those 

15 skilled in the art). It was assumed that no invaded region was present prior to solvent 
injection and that therefore 6 was initially zero throughout the model. Formation of 
the invaded region was triggered by assuming that solvent went exclusively into the 
invaded region at the injection face of the core. After the w h w rh and S w are calculated 
from the above integration, 6 was updated with Eq. (5), and the integration proceeded 

20 to the next time step. The pressure distribution at each time step was then determined 
by integrating Eq. (9) with respect to £ 

Mass Transfer Function 

It was assumed that, as a first approximation, the rate of inter-region transfer 
25 was proportional to the difference between the component's volume fraction in the 
resident and invaded regions: 
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A J = K A X rj - X ij) 



(10) 



where Kj was the mass transfer coefficient for component j [units: time -1 ], and x rj and 
xy = (Sgyj + Spcj)/(1 -S w ) were the volume fractions of component j in the resident and 
invaded regions respectively. In equation (10), the volume fraction difference was the 
driving force for mass transfer and the mass transfer coefficient characterized the 
resistance to mass transfer. With this assumption, equations (6) through (8) became: 



dw rX _ d 



as 



i) 



-Da,(x rl -x n ) 



(11) 



dw 



r2 



dr d% 



4>c 



dw 



r3 



dz 



d_ 

at 



-Da 2 (x r2 -x a ) 



(12) 



-Da 3 (x ri -x a ) 



(13) 



where Daj s ?q$L/u, known as the Damkohler number, was the dimensionless mass 
transfer coefficient. The magnitude of the Damkohler number represented the rate of 
mixing of the component between the invaded and resident regions relative to the 
residence time of fluid in the core. A Damkohler number of zero for all components 
implies no mixing, and high Damkohler numbers implies rapid mixing. 

This model was consistent with the assumption that mixing causes transfer of a 
component from regions of higher concentration to regions of lower concentration, 
thus tending to equalize concentrations between the two regions. 

The mass transfer coefficients may be functions of the local degree of 
miscibility, gridcell geometry, invaded fraction (6), mobility ratio (w), velocity (w), 
heterogeneity, and water saturation (S w ) within the gridcell: 
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Kj = k j {degree of miscibility, gridblock geometry, 0, m, u, heterogeneity, S w ) (14) 

The specific functional dependencies depend on the processes by which the invaded 
and resident fluids mix. Gardner, J. W., and Ypma, J. G. J., "An Investigation of 
Phase-Behavior/Macroscopic Bypassing Interaction in C0 2 Flooding," Society of 
Petroleum Engineering Journal, pages 508-520, October 1984, disclose the effects of 
macroscopic bypassing on mixing in multiple-contact miscible displacement 
processes. The inventors have observed that data presented by Gardner and Ypma 
imply that mass transfer coefficients should be inversely proportional to the time 
required to eliminate subgrid fingers by transverse dispersion: 



(15) 



1 5 where d is the transverse width of the gridcell, D T j is the transverse dispersion 

coefficient of component j, F$ is a parameter accounting for effects of invaded fraction 
and heterogeneity, and C/ y is a constant that may depend on component j. 

As a first approximation, the transverse dispersion coefficient includes 
contributions from molecular diffusion, convective dispersion, and capillary 

20 dispersion. The mass transfer coefficient model incorporates these contributions and 
can be written in dimensionless form as: 



Da } = C l} F 6 



u 



C 2 D oj 4>L + a T (d)<f>L 



d 2 



1 + C y (^—) 



25 



30 



= Da 



Mj 



l + C r (-£-) 

Y max 



(16) 



where D OJ is the molecular diffusion coefficient for component j, aj{d) is transverse 
dispersivity, y max is the maximum gas-oil interfacial tension for immiscible 
displacement, Da M j is the Damkohler number for first-contact miscible displacement, 
and C2 and C r are adjustable constants. The terms in the first bracket are the 
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dimensionless rates of mass transfer due to molecular diffusion and convective 
dispersion, respectively. Molecular diffusion dominates at low velocity and small 
system width, and convective dispersion dominates at high velocity and large system 
width {a T (d) is an increasing function of d). The terms in the second bracket account 



and Da M j are synonymous.) It was assumed for initial testing purposes that the mass 
transfer coefficients were unaffected by mobility ratio and water saturation. 

In multiple-contact miscible and near-miscible displacements, interfacial 
tension depended on the location of the gridcell composition within the two-phase 

1 0 region of the phase diagram; the closer the composition was to the critical point, the 
lower would be interfacial tension. Within the context of the present model, where 
interfacial tension was a measure of the degree of miscibility between solvent and oil, 
the interfacial tension in Eq. (16) was the tension that would exist between vapor and 
liquid if the entire contents of the gridcell was at equilibrium. The following parachor 

15 equation was used to calculate interfacial tension: 



where Pj is the parachor parameter for component j, Xj and yj are the mole fractions of 
component j in the invaded liquid and invaded vapor phases, respectively, Q and £ v are 

20 molar densities of the liquid and vapor and n is an exponent in the range 3.67 to 4. 

A key feature of the mechanistic mass transfer model used in this example was 
that the degree of miscibility between solvent and oil had a significant impact on the 
rate of mixing between the invaded and resident regions. It has been proposed in the 
prior art that immiscible dispersion coefficients of fluids in porous media can be about 

25 an order of magnitude greater than miscible dispersion coefficients under equivalent 
experimental conditions. Therefore, mixing should be more rapid under immiscible 
conditions than under miscible conditions. In the model used in the example, this 
observation was incorporated by including an interfacial tension dependence in the 
calculation of the transverse dispersion coefficient. Since the interfacial tension 



5 



for capillary dispersion. (Note that when C r is zero, i.e., the fluids are miscible, Dctj 




(17) 
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depends on phase behavior through the parachor equation, Eq. (17), the relevant 
parameter in the context of the model was the interfacial tension constant, C r . 

The mass transfer model introduced a number of parameters (e.g., diffusion 
coefficients, dispersivity, interfacial tension) into the predictive model of this 
5 invention that have no counterparts in the Todd-Longstaff mixing model. While these 
additional parameters increase computational complexity, in contrast to the Todd- 
Longstaff mixing model, all parameters of the present inventive model have a physical 
significance that can either be measured or estimated in a relatively unambiguous 
manner. 

10 

Effective Medium Mobility Function 

Percolation theory and the effective medium approximation are known 
techniques for describing critical phenomena, conductance, diffusion and flow in 
disordered heterogeneous systems (see for example, Kirkpatrick, S., "Classical 

1 5 Transport in Disordered Media: Scaling and Effective-Medium Theories," Phys. Rev. 
Lett., 27 (1971); Mohanty, K. K., Ottino, J. M. and Davis, H. T., "Reaction and 
transport in disordered composite media: introduction of percolation concepts," Chem. 
Engng. Sci., 1982, 37, 905-924; and Sahimi, M., Hughes, B. D., Scriven, L. E. and 
Davis, H. T., "Stochastic transport in disordered systems," J. Chem. Phys., 1983, 78, 

20 6849-6864). In the context of flow problems in heterogeneous systems, the effective 
medium approximation represented transport in a random heterogeneous medium by 
transport in an equivalent (effective) homogeneous medium. The inventors have 
observed that the agreement between the effective medium approximation and 
theoretical results is quite good when far away from the percolation threshold. 

25 An effective medium mobility model was generated to evaluate mobilities of 

fluids in a heterogeneous medium. This was done by assuming that the distribution of 
solvent and oil within a region of a gridcell could be represented by a random 
intermingled network of the two fluids. The following analytical expressions for 
nonaqueous phase mobilities were derived by assuming the network to be isotropic 

30 and uncorrected: 
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15 
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(25) 
(26) 



(27) 



The coordination number, z, is a measure of the "branchiness" of the 
25 intermingled fluid networks. Increasing z leads to more segregation of oil and solvent, 
so that solvent breakthrough is hastened and oil production is delayed. The relative 
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permeabilities were evaluated using the saturation of the fluid within its region. The 
effective medium mobility model provided approximate analytical expressions for 
phase mobilities that take into account the relevant properties (invaded fraction, 
heterogeneity, mobility ratio) in a physically sound manner. Results presented below 
5 show that the effective medium mobility model accurately captured the recovery 
profiles in miscible displacements. 

Phase Behavior Function 

A simplified pseudo-ternary phase behavior model was used in the examples of 

10 this invention for the one-dimensional simulator. In this model, the compositions of 
mixtures of solvent and oil were characterized in terms of three pseudocomponents: 
CO2, a light oil component, and a heavy oil component. The two-phase envelope in this 
phase model was described by a quadratic equation, the constants of which were 
determined by the compositions for the plait point and the two termini of the envelope 

15 at the boundaries. While only approximately representing a real system, this phase 
model successfully simulated phase behaviors corresponding to differing degrees of 
miscibility such as first-contact miscible (FCM), multiple-contact miscible (MCM) and 
near-miscible (NM). 

Parameters defining the two-phase envelope used in Examples 1-3 are 

20 summarized in Table 1 . The parameters in Table 1 for the MCM case defined a 

pseudo-ternary phase description of the CC>2-Means crude system at 2000 psia (13,790 
kPa) and 100 °F (37.78 °C). The parameters in Table 1 for the FCM and NM cases 
defined a pseudo-ternary phase description that might be obtained at 100 °F (37.78 
°C) and pressures higher and lower than 2000 psia (13,790 kPa), respectively. The 

25 resident oil composition was predominantly heavy, corresponding to a heavy oil 
fraction of 0.8434 and a light oil fraction of 0.1566. 
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Table 1 



Parameter 


Value 


Vic 


0.99 


V 2G (1-V ia ) 


0.01 


V 3G 


0 


V, L 


0.19197 


V 2L (I-V20) 


0.80803 


v 3L 


0 


rFCM 
V 3P J MCM 
L NM 


0.00 


0.09 


0.36 


rFCM 
V 2P J MCM 
L NM 


0.6372 


0.5472 


0.3072 


v IP 


0.3628 



5 

Referring to Table 1 , the subscripts 1 , 2 and 3 denote solvent, the heavy oil and 
light oil, respectively. Vig and Vil represent the termini of a two-phase envelope. 
Vig and Vil represent the solvent volume fractions in gas and liquid phases 
respectively for the solvent-heavy end mixture. V JP and V3P represent the solvent and 

1 0 light end volume fractions at the plait point. 

Parameters defining the two-phase envelope used in Example 4 (discussed in 
more detail below) are summarized in Table 2. Parameters used in Example 4 defined 
a pseudo-ternary phase description of the CC>2-Wasson crude system at 2000 psia 
(13,790 kPa) and 100 °F (37.78 °C). The data were obtained from Gardner, J. W., 

15 Orr, F. M., and Patel, P. D., "The Effect of Phase Behavior on C0 2 Flood 

Displacement Efficiency," Journal of Petroleum Technology, Nov. 1981, pages 2067- 
2081. The crude oil composition corresponded to a heavy oil volume fraction of 0.72 
and a light oil volume fraction of 0.28. 
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Table 2 



Parameter 


Value 


V,o 


0.97 


v 2G 


0.03 


v 3G 


0 


V, L 


0.23 


V 2L 


0.77 


v 3L 


0 


v 3P 


0.17 


v 2P 


0.48 


Vip 


0.35 



5 Simulation Results 

The input data used in the four example simulations assumed oil-brine relative 
permeability and capillary pressure data characteristic of San Andres carbonate rock. 
Core properties were length = 1 ft (0.3048 m), porosity = 0.19 %, and permeability = 
160 md (0.1579 /mi 2 ). 
10 Example 1 

The coordination number, z, in the effective medium approximation to the 
percolation theory denotes the "branchiness" or connectivity of the network. In the 
context of this invention, z represented finger structure in a gridcell and incorporates 
the effects of properties such as oil/solvent mobility ratio, reservoir heterogeneity, and 

1 5 rock type. In a general way, z may be analogized to the mixing parameter co in the 
Todd-Longstaff mixing model. Fig. 5A shows that increasing z results in reduced oil 
recovery and Fig. 5B sho\ys that increasing z results in earlier solvent breakthrough. 
Both the oil recovery and solvent breakthrough curves are sensitive to the value of z. 
In particular, varying z between two and five reduces oil recovery at 1.5 pore volumes 

20 produced from 93% to 52% and reduces the point at which the produced fluid reaches 
a concentration of 50% solvent from 0.55 to 0.24 pore volumes produced. The MCM 
phase behavior description in Table 1 was used in this example and the Damkohler 
numbers were assumed to be Da\ = 0, Da2 = 0.1, and Da^ = 0.1. The simulation of 
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this example started at a waterflood residual oil saturation of 0.35 and used 25 
gridcells in the one-dimensional model. 

An increase in the value of z in effective medium model produced an effect 
similar to a decrease in the value of the mixing parameter co in the Todd-Longstaff 
5 mixing model; both resulted in increased bypassing of oil (lower recovery) and earlier 
solvent breakthrough. The coordination number z can be assigned values greater than 
or equal to two in the practice of the method of this invention, z = 2 represents flow 
of oil and solvent in series and characterizes a piston-like displacement with no 
fingering or channeling, z —> oo represents flow of oil and solvent in parallel and 
1 0 characterizes a displacement with extensive fingering or channeling. Based on these 
results, z can be expected to be important parameter in matching solvent breakthrough 
and oil production history. 

Example 2 

The Damkohler numbers represent the rate of mixing of components between 
15 invaded and resident regions. Results shown in Figs. 6A through D demonstrate that 
this invention successfully reproduces the correct limiting behaviors. The MCM 
phase behavior description in Table 1 was used in this example and the Damkohler 
numbers were assumed to be Da\ = 0 for the solvent component and Da 2 = Da 3 for the 
oil components. The simulation of this example started at a waterflood residual oil 
20 saturation of 0.35 and used 25 gridcells in the one-dimensional model. 

Fig. 6A shows that when there is no mixing (oil Damkohler numbers = 0), 
the model correctly predicts that there is pure displacement of the oil with no 
exchange of components between regions. In Fig. 6A, curve 30 is the fraction of light 
oil component recovered and curve 31 (which has exactly the same shape as curve 30) 
25 is the fraction of heavy oil component recovered. The light and heavy component 

recovery curves 30 and 31 are identical, which indicates that the composition of the oil 
did not change. 

When there is rapid mixing (oil Damkohler numbers greater than about 5), the 
two regions quickly attain nearly identical composition. Therefore, the results of the 
30 simulation shown in Fig. 6D are effectively identical to those of a conventional single- 
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region model. In Fig. 6D, curve 60 is the fraction of light oil component recovered 
and curve 61 is the fraction of heavy oil component recovered. 

The results shown in Fig. 6D also indicate that as the Damkohler number increases in 

a MCM recovery process, there is increasing fractionation of the light oil component 
5 into the gas phase. Consequently, the light component was preferentially recovered as 

the invading (high-mobility) solvent swept it out and left behind residual oil enriched 

in the heavy component. 

Figs. 6B and 6C show results for intermediate rates of mixing. In Fig. 6B, 

curve 40 is the fraction of light oil component recovered and curve 41 is the fraction 
10 of heavy oil component recovered. In Fig. 6C, curve 50 is the fraction of light oil 

component recovered and curve 51 is the fraction of heavy oil component recovered. 

These figures show that the amount and composition of the oil recovered depends 

strongly on the Damkohler numbers. Thus, the timing of each component's recovery 

could be matched by adjusting the Damkohler numbers. Small changes in oil recovery 
1 5 and matching the produced oil and gas compositions could be accomplished through 

variation of the Damkohler numbers. 

Example 3 

Fig. 7 shows experimental data presented in a paper by Blackwell, R. J., 
Rayne, J. R., and Terry, W. M, "Factors Influencing the Efficiency of Miscible 

20 Displacement," Petroleum Transactions, AIME (1959) 216, 1-8 (referred to 

hereinafter as "Blackwell et al ") for a first-contact miscible flood at different values of 
initial oil/solvent viscosity ratio. The experimental data, which appear as points in Fig. 
7, were obtained using homogeneous sand packs and fluids of equal density (to 
minimize gravity segregation). Experiments were conducted at viscosity ratios of 5, 

25 86, 150 and 375. No water was present in the experiments. 

Also plotted in Fig. 7 are lines that correspond to oil recoveries obtained from 
simulations using the method of this invention in which the initial oil/solvent viscosity 
ratio was set at the experimental value, and the coordination number was adjusted to 
obtain the best possible fit with the experimental data. The Damkohler number was 

30 estimated to be on the order of 1 0" 4 (based on D T = 0.0045 ft 2 /day (4.2 cm 2 /day), $ = 
0.4, L = 6 ft (1.83 m), d= 2 ft (0.61 m), and u = 40 ft/day (12.2 m/day)) and was 
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therefore assumed to be effectively zero. There is thus only one adjustable parameter 
used in the simulation — the coordination number, z. Twenty-five gridcells were used 
in the one-dimensional model. 

Fig. 7 shows excellent agreement between the experimental data of Blackwell 
5 et al and results generated using the method of this invention. In particular, the 
method of this invention successfully predicted the leveling-off of the oil recovery 
after initial breakthrough. Moreover, the agreement with the data points for the 
adverse viscosity ratio displacements was exceptionally good. Since the system 
employed by Blackwell et al. was first contact miscible and dispersion was negligible, 

10 neither phase behavior nor mass transfer played a role in the change in simulated 

recoveries. The agreement with experiment in this instance is therefore a validation 
only of the effective medium model of this invention. 

While the procedure adopted above may be equated with history matching 
field data, for the method of this invention to have predictive capability, it would be 

1 5 necessary to be able to predict the value of z a priori. The choice of z would be 
influenced by the mobility ratio, the reservoir heterogeneity and rock type. Fig. 8 
shows a plot of the z values that were used to obtain the fits with experimental data in 
Fig. 7 as a function of oil/solvent viscosity ratio. As illustrated in Fig. 8, z shows a 
monotonic variation with viscosity ratio. 

20 The results presented in Examples 1 and 3 indicate that the coordination 

number, z, is a key parameter in the practice of this invention since it can be used in 
matching solvent breakthrough and oil production history. Example 2 indicates that 
fine tuning of oil recovery as well as matching the produced oil and gas compositions 
can be accomplished through the mass transfer model. 

25 Using the coordination number, z, and the Damkohler numbers as adjustable 

parameters, and the appropriate phase model for the system under study, the predictive 
model of this invention could be used to match the essential features (including oil 
recovery, injected fluid breakthrough, and produced fluid compositions) of any gas 
injection process. 
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Example 3 indicates that the effective medium mobility model used in the 
method of this invention can be used to describe the fingering and bypassing that is 
prevalent in miscible displacement processes. 
Example 4 

5 Example 4 is presented to demonstrate the utility of the phase behavior and 

mass transfer models. Experimental data presented in papers by Gardner, J. W. 5 Orr, 
F. M., and Patel, P. D. 5 "The Effect of Phase Behavior on C0 2 Flood Displacement 
Efficiency " Journal of Petroleum Technology, pages 2067-2081, Nov. 1981 (referred 
to hereinafter as "Gardner et a/. M ) and Gardner, J. W., and Ypma, J. G. J., "An 

10 Investigation of Phase-Behavior/Macroscopic Bypassing Interaction in C0 2 

Flooding," Society of Petroleum Engineers Journal, pages 508-520, October 1984, 
disclosed the relationship between phase behavior and displacement efficiency (oil 
recovery) for miscible gas injection processes. These papers presented results of 
coreflood experiments on two systems: (i) displacement of Soltrol by CO2 in a first 

15 contact miscible (FCM) system, and (ii) displacement of Wasson crude by CO2 in a 
multiple-contact miscible (MCM) system. Soltrol is a product manufactured by 
Phillips Petroleum Company and Wasson crude is from the Wasson field in west 
Texas. The oil/solvent viscosity ratio was 16 for the CCVSoltrol system and 21 for 
the CCVWasson crude system - close enough so as to make phase behavior the only 

20 major distinction between the two systems. Therefore, for all practical purposes, the 
only reason for any difference in recoveries for the two systems could be attributed to 
the change in phase behavior and macroscopic bypassing (as a result of the changed 
phase behavior). 

Fig. 9 shows the experimental recovery curves obtained for the CCVSoltrol 
25 (curve 70) and C0 2 /Wasson (curve 71) crude systems. The different sets of symbols 
denote data obtained in duplicate coreflood experiments under similar conditions. All 
tests were done in the same Berea core. Ultimate oil recovery efficiency was lower 
for the CO2/ Wasson crude system, as was the rate of recovery. 

Viscous fingering was almost entirely responsible for the shape of the FCM 
30 CCVSoltrol recovery curve 70 while both viscous fingering and phase behavior were 
responsible for the shape of the MCM CCVWasson crude recovery curve 71. To test 
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the influence of fingering on recovery, one-dimensional simulations were first run 
using a conventional single-region model. For the simulations of this example, 
simulation parameters were set to closely match the CC^/Soltrol and COi/Wasson 
crude experimental systems. The C0 2 viscosity was set at 0.063 cp (0.000063 Pa/sec) 
5 in line with data provided by Gardner et ah Soltrol has a normal boiling point range 
equivalent to that of C11-C14, which corresponds to a viscosity of approximately 1.2 cp 
(0.0012 Pa/sec). However, in order to exactly match the experimental oil/solvent 
viscosity ratio of 16, the Soltrol viscosity was assumed to be 1.01 cp (0.00101 Pa/sec). 
Phase viscosities were calculated by the quarter-power blending rule, which is well 

10 known to persons of ordinary skill in the art. 

Experimental gas/oil relative permeability ratios were used in establishing the 
relative permeability-saturation relationship in the simulation. The simulations were 
run with 30 gridcells. The number of gridcells was chosen so as to approximate the 
level of longitudinal dispersion in the experimental systems. In the case of the 

1 5 CCVWasson crude simulation, the phase model was chosen to be the same as the 
experimental one, shown in Table 2. Fig. 9 shows the recovery curves 72 and 73 
obtained from the single-region model simulations along with the experimental data 
(curves 70 and 71). Curve 72 illustrates simulation results of the CC^/Soltrol system 
and curve 73 illustrates simulation results of the CCVWasson system. It is clear from 

20 Fig. 9 that viscous fingering suppresses the rate of recovery of oil. It is also apparent 
that the single-region model provides an inadequate description (both qualitatively as 
well as quantitatively) of the oil recovery in the CCVSoltrol and CC^/Wasson crude 
systems. However, the single-region simulations are in good agreement with slim 
tube experiments (Gardner et al.) in which the effects of bypassing are suppressed. 

25 To evaluate the ability of the method of this invention to simulate the 

experimental coreflood data, the method of this invention was first applied to the FCM 
CC>2/Soltrol system. The parameters z, Da so i V ent, DaMheavy and Da M ii g ht were adjusted so 
as to obtain the best possible fit with the experimental data. DaMheavy was assumed to 
be equal to DaMtight for simplicity. The best fit was obtained for the selection z — 4.5, 

30 Da so i vent = 0, DaMheavy, light = 0.5. Using the same parameters and assuming C r = 10, a 
simulation was carried out using the method of this invention for the CCVWasson 

I:\UR<^-LAW\PATENTS\RES\97095\PROVJSIONAL APPLlCATION.DOC 




-36- 

crude system. All simulation parameters (phase behavior, relative permeability- 
saturation relationship and dispersion level) were set to match the experimentally 
determined values (data obtained from Gardner et al). The viscosity of the oil in the 
simulation was changed to mimic the Wasson crude and an oil/solvent viscosity ratio 
5 of 21 . These results are plotted in Fig. 10. 

In Fig. 10, curves 70 and 71 of Fig. 9 are again shown to compare the 
simulation results, curve 74, of the CCVSoltrol system using the two-region model of 
this invention and simulation results, curve 75 of the CCVWasson crude system using 
the two-region model used in the method of this invention. 

10 The method of this invention did an excellent job of matching the MCM 

CC>2/Wasson using the same parameters that were applied to the FCM CC>2/Soltrol 
crude system. The rationale for keeping z fixed from the CC^/Soltrol simulation is 
that, since the Soltrol and Wasson crude experiments were conducted on the same 
cores (same degree of heterogeneity and rock type), and at virtually the same 

15 oil/solvent viscosity ratio (same mobility ratio), the value of z must remain essentially 
unchanged. Mass transfer coefficients increased from the values used for the best fit 
of the CCVSoltrol system. Physically, this translates into an increase in mass transfer 
rates with reduction in miscibility (FCM to MCM) - as miscibility decreases, capillary 
dispersion increases resulting in higher rates of mass transfer. 

20 In the simulations presented in the foregoing examples, it was assumed that the 

resident region remained a single-phase liquid. However, the composition of the 
resident region may enter into the multiphase envelope if solvent components are 
allowed to transfer into that region, which could be performed by persons skilled in 
the art. This would necessitate an additional flash calculation for the resident region 

25 and the need to specify both vapor and liquid phase permeabilities for that region. 
The Partitioned Node Model used in the method of this invention is 
particularly attractive for use in modeling solvent-flooded reservoirs because all the 
parameters used in the model have a physical significance that can either be measured 
or estimated by those skilled in the art. 
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The coordination number, z, in the effective-medium model can be adjusted to 
match the timing of injected fluid production. It has been observed that z increases 
with increasing initial oil/solvent mobility ratio. 

The constants, C/y, in the mass transfer function can be adjusted to match 
5 individual component production histories. Molecular diffusion coefficients, Z) oy , can 
be estimated with standard correlations known to those skilled in the art. Dispersivity, 
or, and the diffusion constant, C?, will depend on rock properties, and will determine 
scaling from laboratory to field. In most applications, the interfacial tension 
parameter, should be a constant, to good approximation. 

10 The effect of gravity on relative mobilities, which was not addressed in 

foregoing examples, can be also be taken into account by those skilled in the art. For 
example, it may be expected that within a gridcell, the low-density phase would tend 
to segregate to the top of the gridcell and would have a higher effective mobility in the 
upward direction. Anisotropy in permeability was also not considered in the example 

15 simulations. In a 3-D simulation, absence of such anisotropy may tend to overestimate 
flow in the vertical direction. An anisotropic formulation of the effective medium 
model can be incorporated into the model by those skilled in the art, but this would 
significantly increase the complexity of the computations. 

A still another factor that was not considered in the present examples was the 

20 presence of water in the gridcells. In simulating water-altemating-gas (WAG) 

injection, gas would be injected only into the invaded region and water would only be 
injected into the resident region. In this way, formation of the invaded region would 
be triggered only by injection of the high-mobility gas and not by injection of water. 
Water saturation could also have an effect on the oil/gas mass transfer coefficients - 

25 which would typically be incorporated into the model. A transfer function can be 

developed for water by those skilled in the art, so that water can also partition between 
the invaded and resident regions. 

The principle of the invention and the best mode contemplated for applying 
that principle have been described. It will be apparent to those skilled in the art that 

30 various changes may be made to the embodiments described above without departing 
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from the spirit and scope of this invention as defined in the following claims. It is, 
therefore, to be understood that this invention is not limited to the specific details 
shown and described. 
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Symbols 

C/j constant used in describing mass transfer coefficient of component j 

C2 ratio of apparent diffusion coefficient in porous medium to 

molecular diffusion coefficient 
C r interfacial tension (IFT) parameter 

D width of gridcell 

Dctheavy Damkohler number of heavy oil component 

Dcij Damkohler number of component j (includes interfacial tension 

effects) 

Dau g ht Damkohler number of light oil component 

Dqmj Damkohler number of component j for first-contact miscible 

displacement (excludes interfacial tension effects) 
Da so i ven t Damkohler number of solvent 

D 0 j molecular diffusion coefficient for component j 

Dry transverse dispersion coefficient of component j 

FCM First-Contact Miscible 

Fq parameter accounting for effects of invaded fraction and 

heterogeneity 
K permeability 
L core/gridcell length 

M mobility ratio 

MCM Multiple-Contact Miscible 

NM Near-Miscible 
P pressure 
p c capillary pressure 

Pj parachor parameter for component j 

Q volumetric injection rate 

S gs St vapor and liquid saturations in the invaded region 

S w water saturation 

T time 
U velocity 

ViGf Vji pseudo-ternary phase description parameters: solvent volume 

fractions in gas and liquid phases for the solvent-heavy end mixture 

Vjp pseudo-ternary phase description parameter: solvent volume fraction 

at the plait point 

V 3P pseudo-ternary phase description parameter: light end volume 

fraction at the plait point 
V p pore volume 

wj,W2, w 3 volume fraction of the solvent, the heavy fraction of the oil and the 

light fraction of the oil 
Wii, Wi2, Wr3 volume fraction of the solvent and heavy fraction of the oil in the 

invaded region 

Wri, Wr2, w r3 volume fraction of the solvent and heavy fraction of the oil in the 

resident region 
X length 
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Xy 


volume fraction of component j in the nonaqueous portion of the 




invaaea region 


xj, yj 


volume fraction of component j in the liquid and vapor portions of 




trie invaded recxinn 

L11W 111 V ClUv/VJ IVglUll 




volume fraction of comnntient i in the nnnamipniK nortion of* the 

vullUllv llCivllV/11 V/l vvlllUUllvlll 1 111 Lllw 1 1U11CIU UviU L/V-/1 L1V_/11 Lilt/ 




resident region 


z 


coordination number 


OCT 


transverse dispersivity 


R 

H 


dimensionless permeability, = k/uL 


Y 
f 


interfacial tension 


V 

[max 


maximum eas-oil interfacial tension for immiscible displacement 


h 


dimensionless length = x/T, 




molar densities of the HmiiH and vanor 


& <P 


uui vjoil y 


~r • js-j 

'-■J ^ 


mass transfer poeffifient nf rnmnnnpnt i 

IIICIOO 11 CUlJl^/l vUvlllvlblll Ul C'UlllL/v/llt'llt I 




rate of transfer f volume/time^ of component i from the resident to 




the invaded region 


j5 ^iW) ^ile> X roe 


effective mobilities of the vapor phase in the invaded region, the 




liquid phase in the invaded region, and the resident fluid. 




total effective mobility, = X ive + Xn e + Z roe + A w 




mobility of water 




invaded fraction of gridcell 


!S r 


dimensionless time, = ut/<f)L 
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